/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/device/host/operation.hxx"


namespace linal {
  namespace host {
    namespace internal {
      /*!
        Apply FLAME style pivot vector to rows : Left + No transpose
      */
      int apply_pivots_ln(FLA_Obj &p, FLA_Obj &A) {
        if (!A.n) return true;

        FLA_Obj AT, AB,   A0, a1,   A2;
        FLA_Obj pt, pb,   p0, rho1, p2;

        FLA_Part_2x1( A, &AT,
                      &AB, 0, FLA_TOP );
        FLA_Part_2x1( p, &pt,
                      &pb, 0, FLA_TOP );

        while ( FLA_Obj_length( pt ) < FLA_Obj_length( p ) ) {
          FLA_Repart_2x1_to_3x1( AT, &A0,
                                 &a1,
                                 AB, &A2, 1, FLA_BOTTOM );

          FLA_Repart_2x1_to_3x1( pt, &p0,
                                 &rho1,
                                 pb, &p2, 1, FLA_BOTTOM );
          //------------------------------------------------------------
          int *idx = (int*)(FLA_Obj_buffer_at_view( rho1 ));

          if (*idx) {
            Flat_ A_2(A2), a_2;

            // * extract row vector with distance "*idx" from a1
            A_2.extract( a_2, 1, A_2.get_n(), *idx-1, 0 );
            FLA_Swap( a1, ~a_2 );
          }

          //------------------------------------------------------------
          FLA_Cont_with_3x1_to_2x1( &AT, A0,
                                    a1,
                                    &AB, A2, FLA_TOP );
          FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                    rho1,
                                    &pb, p2, FLA_TOP );
        }
        return true;
      }

      /*!
        Apply FLAME style pivot vector to rows : Left + Transpose
      */
      int apply_pivots_lt(FLA_Obj &p, FLA_Obj &A) {
        if (!A.n) return true;

        FLA_Obj AT, AB,   A0, a1,   A2;
        FLA_Obj pt, pb,   p0, rho1, p2;

        // ** range should be different here because the pivot loop
        //    goes from bottom to top
        FLA_Part_2x1( A, &AT,
                      &AB, 
                      FLA_Obj_length( A ) - FLA_Obj_length( p ), 
                      FLA_BOTTOM );
        FLA_Part_2x1( p, &pt,
                      &pb, 0, FLA_BOTTOM );

        while ( FLA_Obj_length( pb ) < FLA_Obj_length( p ) ) {
          FLA_Repart_2x1_to_3x1( AT, &A0,
                                 &a1,
                                 AB, &A2, 1, FLA_TOP );

          FLA_Repart_2x1_to_3x1( pt, &p0,
                                 &rho1,
                                 pb, &p2, 1, FLA_TOP );
          //------------------------------------------------------------
          int *idx = (int*)(FLA_Obj_buffer_at_view( rho1 ));

          // index will be different. p and A range is different.
          if (*idx) {
            Flat_ A_2(A2), a_2;

            // * extract row vector with distance "*idx" from a1
            A_2.extract( a_2, 1, A_2.get_n(), *idx-1, 0 );
            FLA_Swap( a1, ~a_2 );
          }

          //------------------------------------------------------------
          FLA_Cont_with_3x1_to_2x1( &AT, A0,
                                    a1,
                                    &AB, A2, FLA_BOTTOM );
          FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                    rho1,
                                    &pb, p2, FLA_BOTTOM );
        }
        return true;
      }

      /*!
        Apply FLAME style pivot vector to columns : Right + No transpose
      */
      int apply_pivots_rn(FLA_Obj &p, FLA_Obj &A) {
        if (!A.m) return true;

        FLA_Obj AL, AR,   A0, a1,   A2;
        FLA_Obj pt, pb,   p0, rho1, p2;

        FLA_Part_1x2( A, &AL, &AR, 
                      FLA_Obj_width( A ) - FLA_Obj_length( p ), 
                      FLA_RIGHT );
        FLA_Part_2x1( p, &pt,
                      &pb, 0, FLA_BOTTOM );

        while ( FLA_Obj_length( pb ) < FLA_Obj_length( p ) ) {
          FLA_Repart_1x2_to_1x3( AL, AR, &A0, &a1, &A2, 1, FLA_LEFT );

          FLA_Repart_2x1_to_3x1( pt, &p0,
                                 &rho1,
                                 pb, &p2, 1, FLA_TOP );
          //------------------------------------------------------------
          int *idx = (int*)(FLA_Obj_buffer_at_view( rho1 ));

          if (*idx) {
            Flat_ A_2(A2), a_2;

            // * extract column vector with distance "*idx" from a1
            A_2.extract( a_2, A_2.get_m(), 1, 0, *idx-1 );
            FLA_Swap( a1, ~a_2 );
          }

          //------------------------------------------------------------
          FLA_Cont_with_1x3_to_1x2( &AL, &AR, A0, a1, A2, FLA_RIGHT );
          FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                    rho1,
                                    &pb, p2, FLA_BOTTOM );
        }
        return true;
      }

      /*!
        Apply FLAME style pivot vector to columns : Right + Transpose
      */
      int apply_pivots_rt(FLA_Obj &p, FLA_Obj &A) {
        if (!A.m) return true;

        FLA_Obj AL, AR,   A0, a1,   A2;
        FLA_Obj pt, pb,   p0, rho1, p2;

        FLA_Part_1x2( A, &AL, &AR, 0, FLA_LEFT );
        FLA_Part_2x1( p, &pt,
                      &pb, 0, FLA_TOP );

        while ( FLA_Obj_length( pt ) < FLA_Obj_length( p ) ) {
          FLA_Repart_1x2_to_1x3( AL, AR, &A0, &a1, &A2, 1, FLA_RIGHT );

          FLA_Repart_2x1_to_3x1( pt, &p0,
                                 &rho1,
                                 pb, &p2, 1, FLA_BOTTOM );
          //------------------------------------------------------------
          int *idx = (int*)(FLA_Obj_buffer_at_view( rho1 ));

          if (*idx) {
            Flat_ A_2(A2), a_2;

            // * extract column vector with distance "*idx" from a1
            A_2.extract( a_2, A_2.get_m(), 1, 0, *idx-1 );
            FLA_Swap( a1, ~a_2 );
          }

          //------------------------------------------------------------
          FLA_Cont_with_1x3_to_1x2( &AL, &AR, A0, a1, A2, FLA_LEFT );
          FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                    rho1,
                                    &pb, p2, FLA_TOP );
        }
        return true;
      }
  
      /*!
        Apply FLAME style pivot vector to matrix
      */
      int apply_pivots( int side, int trans,
                        FLA_Obj &p, FLA_Obj &A ) {
    
        int val = side*1000+trans;

        switch( val ) {
        case 210400: /* left  + no trans */ apply_pivots_ln(p, A);break;
        case 210401: /* left  + trans    */ apply_pivots_lt(p, A);break;
        case 211400: /* right + no trans */ apply_pivots_rn(p, A);break;
        case 211401: /* right + trans    */ apply_pivots_rt(p, A);break;
        default:
          LINAL_ERROR( check_side( side ),        LINAL_ERROR_CHECK_SIDE );
          LINAL_ERROR( check_trans_only( trans ), LINAL_ERROR_CHECK_TRANS );
          break;
        }

        return true;
      }
    }
  }
}
